Island size distributions in submonolayer growth: successful prediction by mean field 
theory with coverage dependent capture numbers 
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We show that mean-field rate equations for submonolayer growth can successfully predict island 
size distributions in the pre-coalescence regime if the full dependence of capture numbers on both 
the island size and the coverage is taken into account. This is demonstrated by extensive Kinetic 
Monte Carlo simulations for a growth kinetics with hit and stick aggregation. A detailed analysis 
of the capture numbers reveals a nonlinear dependence on the island size for small islands. This 
nonlinearity turns out to be crucial for the successful prediction of the island size distribution and 
renders an analytical treatment based on a continuum limit of the mean-field rate equations difficult. 
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The kinetics of submonolayer nuclcation and island 
growth during the initial stage of epitaxial thin film 
growth has been studied intcnsivelyboth experimentally 
and theoretically (for reviews, see [lMsl)- A good under- 
standing of this kinetics assists in tailoring self-organized 
nanostructures and thin film devices for specific needs. 
Mean-field rate equations (MFRE) Q successfully pre- 
dict important features such as the scaling behavior of 
the density of stable islands with respect to the T = D/F 
ratio of the adatom diffusion rate D and incoming flux 
F They seem to fail, however, to predict correctly 

the number densities n s of islands composed of s atoms, 
i.e. the island size distributions (ISD) Q. In this con- 
nection, Ratsch and Venables [l(3| as well as Evans et al. 
[H addressed a still open question: whether the MFRE 
are successful in describing the precise shape of the ISD, 
if the correct dependence of the capture numbers <7 S (0) 
on both s and the coverage were taken into account. 
The answer to this question is not obvious, since the 
MFRE with correct capture numbers a s (<d) still neglect 
(i) many-particle correlation effects [ll[ , (ii) spatial fluc- 
tuations in shapes and capture zones of islands, and (iii) 
coalescence events that, despite rare in the early-stage 
growth, could have a significant influence. 

Various theoretical approaches have been developed 
in the past for obtaining appropriate analytical formu- 
lae or approximate numerical results for the cr s (0) (for 
details, see [H, [llj and references therein). These ap- 
proaches focus on the low-temperature case with criti- 
cal size i = 1, i.e. the case when already dimers can be 
considered as stable (on a time scale, where the ISD in 
the initial growth regime is formed). The roughest ap- 
proach is to neglect the dependence and to use just 
two numbers, <j\ for the adatoms and an average number 
<7 for all stable islands with s > 2, and to fit these num- 
bers to give best agreement with simulated or measured 
data. Alternatively, simulated capture numbers for var- 
ious s at a fixed coverage have been considered [3| 
and used in the analysis of experiments |17l | . As shown 



in Fig. 1, however, neither of these approaches as well as 
a more sophisticated self-consistent treatment 0, HH is 
successful in providing a good description of the ISD as 
obtained from KMC simulations. A first numerical study 
for computing coverage dependent capture numbers has 
been performed in 3, 3] using a level set method. In- 
tegration of the MFRE with the obtained capture num- 
bers gave quantitative agreement with KMC results for 
the island density N, but the statistics was insufficient to 
achieve conclusive answers with respect to the ISD. For 
taking into account the correlation between s and the size 
of capture zone areas, i.e. that larger islands tend to ex- 
hibit larger capture zones, a generalization of the MFRE 
towards an evolution equation for the .joi nt p robability 
of island size and capture area was set up [13l4l5l |. This, 
however, had to be done at the expense of introducing 
additional parameters for considering nuclcation events 
inside the capture zones. 

In this Letter we compute the capture numbers <r s (0) 
as a function of both the island size s and the coverage 
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FIG. 1: Island size distribution obtained from KMC simula- 
tion in comparison with ISDs calculated from an integration 
of the MFRE using three different approximations for the 

<r s (e). 
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O by performing extensive KMC simulations. We show 
that based on these functions the ISD for growth kinet- 
ics with hit and stick aggregation is well predicted by 
the MFRE in the growth regime before coalescence. We 
discuss simplified forms of the capture numbers a s (Q) 
with respect to predicting the ISD, which could render 
an analytical treatment of the problem possible. 

The MFRE for a situation at low temperatures (no 
re-evaporation) with a critical nucleus of size i = 1 and 
consideration of direct impingement of arriving atoms at 
islands are 

—±- =(1 - Q)F - 2Dam\ - Dm a s n s 

S>1 

— 2F 1 K 1 n 1 — K s n s (1) 

— =Dni (a s _in s _i - a s n s ) 

+ Fk s _iti s _i - Fn s n s , s = 2,3, ... (2) 

These equations refer to the growth regime, where coa- 
lescence events of islands should be negligible, and it is 
presumed that only single adatoms are mobile and that 
atom movements between the first and second layer can 
be disregarded. Moreover, adatoms arriving on top of an 
island are not counted, i.e. s in a strict sense refers to 
the number of substrate sites covered by an island (or 
the island area). Accordingly, the deposition flux F of 
adatoms in Eq. ([XJ) has to be restricted to the uncov- 
ered fraction (1 — 9) of the substrate area. The terms 
2Dain\ and FK\n\ describe the nuclcation of dimers due 
to attachment of two adatoms by diffusion and due to di- 
rect impingement, respectively. The term Dn\O s n s de- 
scribes the attachment of adatoms to islands of size s > 1 , 
and FK s n s the direct impingement of deposited atoms to 
boundaries of islands with size s. Dividing Eqs. ([TJ and 
([2]) by F leads to evolution equations with the coverage 
= Ft as independent variable and to a replace of D by 
r = D/F on the right hand side. 

Our KMC simulations are performed with an exact 
continuous-time algorithm and periodic boundary condi- 
tions for "hit and stick" aggregation on a square lattice 
with LxL = 8000 x 8000 sites. The lattice constant is set 
to unity. To calculate the capture numbers a s at the cov- 
erage 0, we use the following procedure: Each simulation 
run is stopped at coverage and the number densities 
n s = N s /L 2 , s = 1,2,... are determined, where N s are 
the numbers of monomers (s = 1) and islands (s > 1). 
Then the simulation is continued for a long time interval 
T without deposition and the following additional rules: 
(i) when an adatom is attaching to an island of size s > 1, 
a counter M s for such attachments is incremented and the 
adatom thereafter repositioned at a randomly selected 
site of the free substrate area (i.e. a site which is neither 
covered nor a nearest neighbor of a covered site); (ii) 
when two adatoms form a dimer, a counter M\ for these 








FIG. 2: (a) Dependence of the capture numbers a s (Q) on s 
for four different fixed coverages; the inset shows the corre- 
sponding k s (0). (b) The coefficients a and b of the asymptote 
cr s (Q) ~ a(Q)s + b(Q), and (c) o\ and a as functions of 6. For 
a convenient extraction of the data in (b) and (c) the follow- 
ing fit function can be used (solid lines): a — 0.103 exp(5.60), 
b = 3.85 - 1.1 exp(7.269), oi = -4.5 + 6.55 exp(3.059), and 
er(0) = -6.8 + 9.8exp(9.39). 

nucleation events is incremented and the two adatoms 
thereafter repositioned randomly as described in (i). In 
this way a stationary state is maintained at the cover- 
age 0. Using the counters, the mean times t s = T/M s , 
s = 1,2,..., for the respective nucleation and attach- 
ment events are determined. Given these times, the cap- 
ture numbers a s are calculated by equating Da s nin s , 
s = 1, 2, . . . with 1/t s , yielding a s = l/[Dn\n s T s \. Aver- 
aging the er s over many simulation runs (configurations) 
finally gives c s (0). The K S (Q) are determined from the 
lengths of the islands boundaries, which arc simultane- 
ously monitored during the simulation and averaged for 
each size s. 

Overall the functions cx s (0) and k s (0) were obtained 
for 57 different values in the range 0.005-0.2 and a 
large number of island sizes for each value of 0, ranging 
up to 1000 values for the largest 0. The typical number 
of nucleation/attachment events for each value was 
10 s . 

Figure shows results for cr s (0) as a function of 
s for four different fixed at V = 10 7 . For large s 
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FIG. 3: (a) Island size distributions for three different cover- 
ages at fixed r = 10 7 , and (b) scaled ISDs for four different 
r. The inset in (b) shows the scaled ISDs in a double-linear 
representation. The symbols mark the results from the KMC 
simulations and the solid lines the results obtained from in- 
tegrating the MFRE equations {TJ and ((2j) with a 3 (Q) and 
k s (0) determined by KMC simulations (see text). 



we find a linear dependence of er s (6) on s at all cov- 
erages, which can be explained [l[ by noting that the 
a s (Q) become proportional to the mean capture zone ar- 
eas A s . Since a double-sized capture zone gives on av- 
erage rise to a double-sized island, it holds A s ~ s and 
hence a s ~ s. The asymptotic behavior can be described 
by fT s (©) ~ a(0)s + 6(6), where the slope a(0) is an 
increasing and the offset 6(0) a decreasing function of 0, 
see Fig. (2b). For small s, a nonlinear dependence of <t s (0) 
on s is found. As shown in the inset of Fig. [5^) , the direct 
capture numbers k s (0) have also a linear dependence on 
s but are approximately independent of 0, i.e. k s (0) ~ 
s. In Fig. [2};) we show the capture number ci(0) re- 
lated to nucleation events and the mean capture number 
ct(0) = YZU<r.(Q)n 3 /N, where N = £~ 2 n s . Thcsc 
functions are important when considering the scaled cap- 
ture numbers a s (Q)/a(Q) as function of the scaled island 
size s/s(Q), where s(0) = J2T=2 sn *l N - 4.7 + 818 at 
T = 10 7 here [H). 

By combining the linear function for large s with a 
polynomial at small s to take into account the nonlincar- 
ity, we fitted the curves in Fig. [3^) and used these fits to 
integrate the MFRE ^ and ©. The data for the rcsult- 
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FIG. 4: Island size distribution for 6 = 0.1 from the KMC 
simulation in comparison with the MFRE results when using 
different simplifications with respect to the functional form of 
the capture numbers (see text). 



ing MFRE-ISDs in Figs. [3^,, b are one of our key findings. 
As shown in Fig. EJa., the MFRE-ISD (solid lines) is for 
all coverages in the precoalescence regime in excellent 
agreement with the corresponding KMC-ISD (symbols) 
obtained from the KMC simulations. A variation of T 
does not affect the quality of agreement, as can be seen 
from Fig. |3Jd, where we plot the scaled ISDs n s s 2 /Q ver- 
sus s/s for a fixed coverage = 0.1 and four different T. 
Moreover, one can infer from this figure that the scaled 
ISDs tends to approach a limiting master curve when 
T — > oo. For comparison with earlier results in the liter- 
ature, we show in the inset of Fig. [SJd the scaled ISDs in 
the more common double-linear plot instead of the linear- 
log representation used otherwise in Figs. [TJ [3k, b, and|U 
We chose this linear-log representation to show that the 
MFRE capture the behavior also correctly in the wings 
at very small (s C s) and very large island sizes (s 3> s). 
In fact, the agreement is seen over about four orders of 
magnitude of n s in Fig. [3Ji. This demonstrates that the 
approximations involved in the MFRE are appropriate to 
predict the ISD with high accuracy for the hit and stick 
aggregation considered here. 

So far we have used the complete functional form for 
cr s (0) and k s (Q). The question arises whether all details 
seen in Fig. [5^) are necessary with respect to a good 
prediction of the ISD. To this end we discuss the following 
simplifications: (i) all k s arc set to zero, (ii) the a s (Q) 
are replaced by eri(0) for s = 1 and er(0) for s > 2 
(and analogously for the k s (0)), and (iii) the asymptotics 
er s (0) ~ a(9)s + 6(0) is used for all s > 2, while we keep 
the cti(0) (again the analogous procedure is used for the 

«.(e)). 

Figure [4] shows the MFRE-ISD resulting from these 
simplifications. Neglecting the k s in Eqs. ([l} and (|2|), 
the ISD is again well predicted, see the dashed line. For 
increasing T the agreement becomes even better (not 
shown). When neglecting the s-depcndence (case (ii)) 
the MFRE-ISD has a maximum still close to the KMC- 
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ISD, but its width is much smaller than that of the KMC- 
ISD. The width of the the respective scaled distribution 
tends to zero for r — > oo. Let us remind that we already 
showed in Fig. [1} that a full neglect of the dependence 
also does not yield a good ISD. In case (iii) the MFRE- 
ISD is also poor in comparison with the KMC-ISD. The 
MFRE-ISD shows a second maximum at s = 2, which is 
caused by the fact that the linear relationship underesti- 
mates the 02(6) value, leading to a higher lifetime and 
correspondingly larger concentration of dimcrs. Gener- 
ally speaking, a linear relationship between (t s (Q) and s 
does not cover the small s behavior but, as one would 
expect, it gives a fair account of the shape of the ISD for 
large s. 

In summary, we have demonstrated that an integration 
of the standard MFRE with coverage-dependent capture 
numbers yields an MFRE-ISD that for hit-and-stick ag- 
gregation is in excellent agreement with the KMC-ISD. 
The full dependence of the capture numbers on both the 
island size and the coverage was determined from exten- 
sive KMC simulations and the functional form was ana- 
lyzed in detail. Despite the fact that a linear dependence 
on the island size holds over almost the entire s-rangc, 
the nonlinear behavior is crucial for a good account of the 
ISD. This implies that it will be difficult to find simple 
functions, which one could use in an analytical contin- 
uum approach for the scaled ISD [l6j |. 

We thank H. Brune, W. Dieterich, and J. Evans for 
very valuable discussions. 



[5] 
[6] 



* Electronic address: mario.einax@tu-ilmenau.de 

* Electronic address: philipp.maass@uni-osnabrueck.de 
URL: |http: //www, st atphys . uni-osnabrueck. de 

[1] J. W. Evans, P. A. Thiel, and M. C. Bartelt, Surf. Sci. 

Rep. 61, 1 (2006). 
[2] T. Michely and J. Krug, Islands, Mound and Atoms 

(Springer, Berlin, 2004). 
[3] H. Brune, Surf. Sci. Rep. 31, 121 (1998). 
[4] J. A. Venables, Phil. Mag. 27, 679 (1973); J. A. Venables, 

G. D. T. Spiller, and M. Hanbiicken, Rep. Prog. Phys. 

47, 399 (1984); J. A. Venables, Surf. Sci. 299, 798 (1994). 



[7; 

[8 
[9 
[10 

[11 

[12 
[13 
[14 

[15 
[16 
[17 

[18 
[19 
[20 



D. Kandel, Phys. Rev. Lett. 78 499 (1997). 

B. Mueller, L. Neldelmann, B. Fischer, H. Brune, and 
K. Kern, Phys. Rev. B 54, 17858 (1996); H. Brune, G. 
S. Bales, J. Jacobsen, C. Borogano, and K. Kern, Phys. 
Rev. B 60, 5991 (1999); J. V. Barth, H. Brune, B. Fis- 
cher, J. Weckesser, and K. Kern, Phys. Rev. Lett. 84, 
1732 (2000). 

M. Einax, Z. Ziehm, W. Dieterich, and P. Maass, Phys. 
Rev. Lett. 99, 016106 (2007); W. Dieterich, M. Einax, 
and P. Maass, Eur. Phys. J. 161, 151 (2008). 
M. Einax, W. Dieterich, and P. Maass, J. Appl. Phys. 
105, 054312 (2009). 

G. S. Bales and D. C. Chrzan, Phys. Rev. B 50, 6057 
(1994) 

C. Ratsch and J. A. Venables, J. Vac. Sci. Technol. A 21, 
S96 (2003). 

When starting from a many-particle master equation for 
the kinetics, nonlinear terms oc nm s , cf. Eq. ([2]), would 
arise from factorizing non-equlibrium correlators of occu- 
pation numbers. 

P. A. Mulheran and J. A. Blackman, Phys. Rev. B 53, 
10261 (1996). 

P. A. Mulheran and D. A. Robbie, Europhys. Lett. 49, 
617 (2000). 

J. G. Amar, M. N. Popescu, and F. Family, Phys. Rev. 

Lett. 86, 3092 (2001); M. N. Popescu, J. G. Amar, and 

F. Family, Phys. Rev. B 64, 205404 (2001). 

J. W. Evans and M. C. Bartelt, Phys. Rev. B 63 235408 

(2001). 

M. C. Bartelt and J. W. Evans, Phys. Rev. B 54, R17359 
(1996). 

M. C. Bartelt, A. K. Schmid, J. W. Evans, and R. Q. 
Hwang, Phys. Rev. Lett., 81, 1901 (1998); M. C. Bartelt, 
C. R. Stoldt, C. J. Jenks, P. A. Thiel, and J. W. Evans, 
Phys. Rev. B 59, 3125 (1999). 

F. G. Gibou, C. Ratsch, M. F. , S. Chen, and R. E. 
Caflisch, Phys. Rev. B 63, 115401 (2001); 
F. Gibou, C. Ratsch, and R. Caflisch, Phys. Rev. B 67, 
155403 (2003). 

From previous investigations by Amar et ai, Phys. Rev. 
B 50, 8781 (1994), a linear relationship of the mean size 
son9 for hit and stick aggregation is known to be valid 
in the aggregation regime (where the number of stable 



islands N 



7 



-1/3 



stays almost constant and s ~ r 1 / 3 © 2 
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